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ABSTRACT 


The design and evaluation of a feedback control system via frequency response methods relies heavily 
upon numerical methods. In application, one can usually develop low-order simulation models 
which for the most part are devoid of numerical problems. However, when complex feedback inter- 
actions, for example, between instrument control systems and their flexible mounting structure, 
must be evaluated, simulation models become moderate to large order and numerical problems 
become common. The fundamental reason for this is that as system order enters the range of about 
20 to 200, many popular tried and true methods are subject to severe numerical problems. An 
attempt is made in this paper to summarize a large body of relevant numerical error analysis litera- 
ture in a language understandable to nonspecialists. The intent is to provide engineers using simula- 
tion models with an “engineering feel” for potential numerical problems without getting intertwined 
in the complexities of the associated mathematical theory. Guidelines are also provided by suggesting 
alternate state-of-the-art methods whicii have good numerical evaluation characteristics. 




CONTENTS 


Page 


ABSTRACT iii 

INTRODUCTION 1 

BACKGROUND 2 

POLYNOMIAL AND RATIONAL FORMS FOR TRANSFER 

FUNCTIONS 4 

NON-OSCILLATORY SYSTEM 6 

DAMPED OSCILLATORY SYSTEM 6 

PROCEDURE 6 

SYSTEM EQUATIONS, PRIMARY DATA 16 

REDUCED ORDER SYSTEM EQUATIONS 23 

CONCLUSION 26 

REFERENCES 27 

APPENDIX A-THE METHOD OF LEVERRIER A-1 


V 




PITFALLS AND GUIDELINES 
FOR THE NUMERICAL EVALUATION 
OF MODERATE-ORDER SYSTEM 
FREQUENCY RESPONSE 


Harold P. Frisch 

Goddard Space Flight Center 
Greenbelt, Maryland 


INTRODUCTION 

The NASA large space system technology (LSST) program is attempting to develop an interactive 
integrated analysis capability (lAC). The objective is to provide engineers with an interdisciplinary 
analysis capability supportive of static, time, and frequency domain design and performance evalua- 
tion methods. In order to achieve this objective, standard analysis methods must be critically 
evaluated to determine their suitability for moderate- and large-order system response analysis. 
Unsuitable methods must be flagged and improved methods developed. 

Most frequency domain analyses are adequately carried out with low- to moderate-order systems; 
i.e., systems of order 20 or less. Unfortunately, low- to moderate-order systems are not always 
adequate; for example, in assessing the cross-coupling effects between two or more on-board in- 
strument control systems via spacecraft system flexibility. Systems such as this require moderate- 
to large-order system models. Within the context of this report, systems of order between about 20 
and 200 are considered to be moderate to large order. In practical application, systems of larger 
order are usually not required if good engineering judgment is used in setting up the simulation 
model. When systems of an order higher than 200 are required and standard algorithms are used, 
resultant computing problems will normally be inordinate due to excessive run-time, memory 
requirements, and round-off error buildup. Solutions in these cases are usually only possible if 
specialized techniques are carefully used. 

This paper is primarily concerned with the numerical methods used to set up and evaluate the 
“transfer function” of single variable feedback control theory and the “transfer function matrix” of 
multivariable feedback control theory. The remarks to be made apply to transfer functions written 
with respect to either the Laplace transform variable s (continuous time systems) or to the Z-trans- 
form variable z (discrete and sampled data systems). 

It will be shown, by presenting the results of a few simple numerical experiments, that commonly 
used numerical methods, while excellent for low-order systems, begin to break down in an unpre- 
dictable manner as system order increases. To avoid this problem, alternate methods are presented 
which can be reliably used for moderate- to large-order systems. 
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It is assumed that transfer functions stem from a physically realizable, well-conditioned system. A 
well-conditioned system is defined as one for which slight perturbations in the system’s physical 
characteristics (mass, geometry, feedback gains, etc.) produce correspondingly slight changes in 
resultant frequency response characteristics. An unstable system is a well-conditioned system if 
slight perturbations in its physical characteristics make it slightly more or less unstable. Ill-condi- 
tioned systems are of minimal importance in application, and therefore are not discussed herein. Ill- 
conditioned systems should be redesigned rather than analyzed. 

The intent of this paper is to summarize, from the theory of condition, information which can be 
used to expose potential pitfalls, and to provide guidelines for selecting numerically reliable methods. 

BACKGROUND 

The historical starting point for the subject of numerical error analysis is usually marked by Refer- 
ence 33, the 1947 von Neumann and Goldstine paper “Numerical Inverting of Matrices of High 
Order.” It is an interesting historical note that in 1947 “high order” referred to matrices of order 
greater than or equal to 10, and “unusually large order” referred to order greater than 100 (Refer- 
ence 33, pages 1022 and 1031). In addition to its historical interest, the discussion in its first 
chapter on the sources of error in a computation provides the framework from which a discussion 
of “condition” can conveniently begin. Four primary sources of error are enumerated: (I) Approxi- 
mations implied by the mathematical model; (2) Errors in observational data, (3) Finitistic approxi- 
mations to transcendental and implicit mathematical formulations; and (4) Errors associated with 
computing instruments (e.g., round-off errors in digital computation). 

In their paper, von Neumann and Goldstine state that their follow-on development pertains only to 
errors of the fourth type. Their objective was to establish rigorous error estimates with respect to 
fbced point digital computation round-off error for the matrix inversion problem. This they do via 
a long and difficult development, which is reviewed in Reference 34 by Wilkinson. It is pointed out 
therein that the complexity of the von Neumann and Goldstine development effectively stymied 
follow-on research in numerical error analysis for the following decade. This situation existed until 
Givens, in Reference 1 3, recognized that round-off errors generated during the course of computa- 
tion could be interpreted as equivalent perturbations on observational data. Furthermore, associated 
perturbation bounds could be established and computed quantities taken as exact solutions to the 
problem defined by the perturbed set of observational data. Nearly all recent analysis makes use of 
this approach. It has the advantage that round-off errors are placed on the same footing as errors in 
input data specification, the effect of which is usually considered during the course of a system 
analysis. This approach is generally referred to as backward error analysis. 

In concept, condition starts with a computing problem. It then asks the following question: If 
the initial data to the computing problem is perturbed by a given amount, how much will the com- 
puted solution differ from the exact solution to the unperturbed problem? If small perturbations in 
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the initial data lead to correspondingly small perturbations in computed results the computing 
problem is considered to be “well-conditioned.” If small perturbations in the initial data lead to 
large perturbations in computed results the computing problem is considered to be “ill-conditioned.” 
Beyond this point, theory becomes complicated by subtle differences in approach and by the 
recognition that for multi-output computing problems; e.g. eigenanalysis, part of the output may be 
well-conditioned, wlule the rest may be either ill-conditioned or almost ill-conditioned. 

A general theory of condition is presented in Reference 25 by Rice. In that reference. Rice at- 
tempts to analyze the way in which an uncertainty (perturbation) in the initial data propagates to 
an uncertainty in the solution of the system of equations. While the results presented by Rice can 
be applied to computational problems, they cannot be construed to form a theory of the condition 
of a computation (see Reference 25, page 288). 

Before one erroneously links round-off error and ill-conditioning together, it is important to point 
out that ill-conditioning exists independently of round-off error effects and that its intrinsic nature 
is unaffected by various schemes for round-off error control. The link between round-off error and 
ill-conditioning comes about after the decision is made to view the accumulated effects of round-off 
error as a perturbation on observational data. 

Wilkinson applies this point of view in his discussions on ill-conditioned problems in Reference 
35. On page 28 he defines a computing problem as ill-conditioned if very small relative perturba- 
tions in the parameters make comparatively large errors in the solutions. 

Steward, in Reference 28, page 76, defines a “stable algorithm” as an algorithm which yields a 
solution which is near the exact solution to the slightly perturbed problem. For ill-conditioned 
problems, the exact solution of the perturbed problem may not be near the exact solution of the 
unperturbed problem and thus, even though a stable algorithm is used to solve an ill-conditioned 
problem, the computed solution and the exact solution need not agree at all. When applied to a 
well-conditioned problem, errors introduced by the stable algorithm are no more than are warranted 
by the data. 

The motivation for this paper stems primarily from the need to obtain transfer functions of moder- 
ate- to large-order systems and the following statement by Wilkinson in Reference 34, page 565. In 
effect it says that certain instinctive methods for setting up and evaluating transfer functions should 
be avoided: 

“Perhaps the most important lesson that has been learned is the rather negative one, that in 
general it is undesirable to reduce implicit polynomial equations, particularly determinantal 
equations, to explicit polynomial form. Such a reduction is superficially attractive both 
because of the formal simplicity of the explicit form and because there is a natural tendency 
to think that ‘we know all about polynomials.’ Commonly a catastrophic worsening of the 
condition of the problem takes place.” 
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The objective of this paper is to explore the dilemma brought about by the fact that the instinctive 
and most commonly used method for obtaining transfer functions is “by reducing implicit polyno- 
mial equations, particularly determinantal equations, to explicit polynomial form.” This exploration 
is to be done via numerical examples designed to provide the engineer with some insight as to how 
far instinctive methods can be reliably used, and then to define what should be done when the 
instinctive methods begin to numerically fail. 

POLYNOMIAL AND RATIONAL FORMS FOR TRANSFER FUNCTIONS 

Virtually all frequency-response methods (see Reference 1 9) require transfer functions to be defined 
as rational functions in the form: 


R(x) = 


Pm 

Q„ « 


( 1 ) 


where the degree m of the numerator polynomial (x) is less than or equal to the degree n of the 
denominator polynomial Q„(x). Multivariable methods require transfer function matrices to be 
defined as rectangular matrices of rational functions. These normally are expressed in the form; 


G (x) = 


N(x) 
d (x) 


( 2 ) 


where N(x) is a rectangular matrix of polynomials and d(x) is the polynomial which is the least 
common denominator of all elements in G(x). 

If the frequency response analysis is pursued analytically, the polynomials contained in equa- 
tions 1 and 2 may be expressed in a variety of equivalent forms, whichever is most convenient for 
the follow-on analytic development. However, if the need for numerical evaluation is on the 
bottom tine, it must be recognized that the evaluation properties of equivalent polynomial forms 
differ drastically relative to each other. (Refer to Rice, Reference 24 and Gautschi, Reference 12 
for an in-depth study.) In summary. Rice, Gautschi, Wilkinson (in Reference 35), and several other 
authors conclude that polynomials expressed in the “root product” form, that is: 


Pn W " ~[J (X-Z.) (3) 

i = 1 
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are comparatively well-conditioned relative to polynomials expressed in the “nested” or “power 
series” form, that is: 


P„ (x) = (• • *((A^ x + A^ j)x + A_j 2 >x+ • • • + Aj)x + Aq 



(4) 


For frequency response analysis, this is not the end (see Wilkinson, Reference 35, page 47). In 
nearly all frequency response problems, neither the zeros of equation (3) nor the coefficients of 
equation (4) can be considered as primary system definition data. Therefore, the problem of 
determining whether or not these parameters have been evaluated with an accuracy sufficient for 
the numerical evaluation of the system’s “true” frequency response characteristics must also be 
considered. 

The following portion of this section explores this problem by presenting the results of a numerical 
experiment. The experiment and data presentation format has been structured to provide insight; 
the actual numeric magnitudes are of minimal importance. In fact, if the experiment were to be 
repeated on a different computer with different precision, a one-for-one correlation would not be 
expected. The important point is that computational problems can develop rapidly, as system order 
increases for certain classes of problems. 

Methods used to compute roots for the root product form will not be explored. If roots are obtained 
via the EISPACK eigenanalysis methods provided in Reference 27 and if primary system definition 
data is used to define associated matrix coefficients, then, since all algorithms given therein are 
stable, errors involved in the computation of eigenvalues are no more than are warranted by the 
system definition data. Refer to the relevant papers collected in Reference 36 for details; nearly all 
EISPACK routines are based upon the ALGOL procedures presented and developed therein. 

The power series form for polynomials is a favorite in controls analysis. Coefficients are easily 
computed and zeros obtainable via application of a variety of stable, well-established polynomial 
root-finding algorithms; e.g. RPOLY (Reference 15). Unfortunately, the use of a stable algorithm is 
no guaranty that accurate results will be obtained for an ill-conditioned computing problem. As 
system order increases, polynomial coefficients can only be computed with a finite degree of 
accuracy on any computing system, i.e., to machine precision. This limit on achievable accuracy, 
when viewed as a perturbation on input data, leads to errors in computed results which far exceed 
that warranted by the imprecision associated with the computed polynomial coefficients. 

In order to illustrate this point in a manner which will provide some physical intuition, two test 
problems are chosen; first, a hypothetical non-oscillatory system; and second, a hypothetical 
damped oscillatory system. 
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Non-Oscillatory System 


Let: 


(X + XjHX + X^H- • -)(X + X^) = 0 


(5) 


where X. = 1, 2, • • N define the first system. 

Damped Oscillatory System 
Let: 


(X^ + 2fcOj X + oj])(X^+ 2?co^ X + (• • •) (X^ + 2foj^ X + cop - 0 (6) 


where co. = 1, 2, * • *, N and 0 < f < 1 define the second. 

In actuality, many values of f have been examined, however results are only presented for f = 
.0001 , .9, and 1 .0. For insight purposes, this is sufficient. Results for other values of f are obtain- 
able, approximately, from an exponential curve fitted between the three relevant points provided. 

Procedure 

Step 1. Utilizing floating point double precision arithmetic, compute the coefficients to the power 
series equivalents of equations 5 and 6. UseN = 2,3, • • *, 25 forequation 5, and N = 2, 3, • ’ *,19 
for equation 6 and the values of f stated above. 

Step 2. Transfer with full precision the coefficients for each polynomial into the zero-finding 
algorithm. The algorithm used was an in-house variation of the popular Jenkins and Traub algo- 
rithm RPOLY given in Reference 15 and discussed in Reference 16. To check root quality, all 
inputted coefficients must be reproducible to 7 digits accuracy with the computed zeros. If not, an 
error flag is automatically set prior to subroutine exit. In all cases presented herein, successful 
quality checks were obtained. 

Step 3. Compute the number of digits in agreement between computed and exact solutions via the 
equation: 


d = 


-log,o 


e - c 


(7) 
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where: 


e = exact solution 

c = computed solution 

d = number of digits in agreement 

All computed zeros are arranged in increasing order of magnitude and compared with the exact 
solutions which are also arranged in increasing order of magnitude. 

All results are presented in Tables 1 through 8, in matrix format. The element d^^j in row n column 
j is to be interpreted as meaning that “the j-th zero in increasing order of magnitude of the polyno- 
mial of degree n has been computed with d^^^ digits of accuracy.” If d^^. > 15.7 this is to be inter- 
preted to imply exact agreement. All results were obtained from an IBM 360-91, using 64 bits 
for double precision computation. 

A full set of analogous results was also obtained by using the EISPACK (Reference 27) eigenvalue 
finding algorithms. In theory, the eigenvalues of the companion matrix (Frobenius form) to the 
associated characteristic polynomial equal the zeros of the polynomial. In actuality, eigenvalues 
computed via the QR-algorithm were inferior in quality to the zeros computed via the polynomial 
root-finding algorithm. However, they were surprisingly good, considering the enormous spread in 
magnitude of the non-zero elements in the companion matrix. This observation is in agreement 
with round-off error analysis theory (see Reference 37, pages 432-434 and Reference 14, pages 
147-148). 

Tables 1 and 2 present results associated with equation 5. This is a favorite example in the literature; it 
has been studied by Gautschi (Reference 12, page 410) and Wilkinson (Reference 37, page 418). In 
Table 1 it is noted that there is a steady deterioration in computational accuracy, as system order 
increases from 2 through 25. Twenty (20) marks the point at which computed results deviate 
significantly enough from the exact for false conclusions to be a real problem. Up to N = 22, all 
computed roots are real, as they should be, as shown in Table 2. For N = 23 and beyond, some 
adjacent real roots are computed incorrectly as complex pairs. This phenomenon is well-known, 
and reported also in Reference 37. 

Tables 3 and 4 present results associated with equation 6 for f = .0001. This is the case of a nearly 
undamped system. It is evident from the results presented that both real and imaginary parts of all 
zeros are computed to a surprisingly high degree of accuracy for polynomials up to degree 38. An 
important point which is not obvious from the tables is that in all cases the real part of each zero is 
computed with the correct sign. 

Tables 5 and 6 present results associated with equation 6 forf = .90. These results are intended to 
provide data at the approximate point at which a breakdown in computational ability becomes 
significant, that is, when system roots approach coincidence in the complex plane or, in engineering 
terms, when the damping associated with system frequencies approaches “critical.” 
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Table 1 - Root 1 , 2, . . N 
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(20) 

z 
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15.7 

15.7 
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15. 1 

15. 7 

15.5 


















4 

1 

15. 1 

1 4. 6 

14.5 

15.0 

















5 

1 

15.7 

14. 9 

14.7 

15.6 

15.2 
















6 

1 

15.2 

13.8 

13 .1 

12.7 

12.7 

13.2 















7 

1 

15.2 

13.3 

12.6 

12.2 

12.1 

12.2 

12.8 














a 

1 

14.4 

13.6 

13.0 

12.5 

12.6 

12.9 

12.5 

1 2.9 













9 

1 

15.4 

14.1 

12.6 

1 1 .d 

1 1.2 

10.9 

10.8 

11.0 

1 1.6 












\ 0 

1 

15. a 

13.7 

12.4 

11.5 

L 0 .9 

10. 6 

10.4 

1 0.5 

10.9 

11.6 











1 1 

1 

15.0 

13.3 

1 1.9 

11. U 

10.4 

10.0 

9.7 

9.6 

9.7 

10. 1 

10.8 










12 

1 

14.V 

12. 7 

11.3 

10.3 

9.6 

9. 1 

8.8 

a.6 

8.7 

8.9 

9*4 

10.3 









13 

1 

15.0 

12.8 

1 1 .4 

10.5 

9.9 

9.6 

9.9 

8.8 

8.4 

8.4 

8.6 

9.0 

9.9 








1 4 

1 

14.0 

12.6 

11.1 

10 .o 

9.1 

d.5 

8. 1 

7.8 

7.8 

8.0 

9.1 

8.6 

8.8 

9.6 







15 

1 

14.9 

12.4 

11.1 

10.0 

9.2 

8.7 

8.4 

9.4 

7.8 

7.4 

7.3 

7.4 

7.8 

6.3 

9.2 






1 6 

1 

14.5 

12.5 

1 1 .0 

9.9 

d.9 

e. 1 

7.5 

7.1 

6.8 

6.7 

6.7 

6.9 

7.2 

7.5 

6*0 

6.6 





1 7 

1 

14. 7 

12.9 

10 .8 

9.4 

8 .4 

7.5 

c . b 

6.3 

5.9 

5*6 

5.4 

5.3 

5.4 

5.7 

6.2 

6*8 

7.9 




1 a 

1 

14.4 

12.5 

10.9 

10.0 

e.7 

7.3 

6.3 

5.6 

5.0 

4.6 

4.3 

4.1 

4.1 

4.2 

4.4 

4.6 

5.5 

6.4 



1 9 

1 

1 4 . 5 

12.2 

10.5 

9.1 

d.o 

7.1 

6.4 

6.9 

5*6 

5.4 

5.2 

4.9 

4.5 

4.3 

4.3 

4.5 

4.8 

5.3 

6.3 


20 

1 

14.4 

12.3 

9.7 

d . 6 

7.7 

7.0 

6.4 

5.8 

5.2 

4.8 

4.5 

4.7 

4.3 

3.7 

3.5 

3.5 

3.6 

4.0 

4.5 

5.5 

21 

1 

13. a 

11.6 

9.4 

d . 1 

6.9 

6.0 

5.1 

4.4 

3.8 

3.2 

2.8 

2.4 

2.2 

2.0 

2.0 

2.0 

2.2 

2.6 

3.0 

3.7 

2 1 

2 1 

4.7 




















22 

1 

14.5 

11.6 

9.3 

0.0 

6.8 

5.9 

5.3 

4.9 

4.7 

3*6 

3.0 

2.7 

2.7 

4.9 

2.3 

1.9 

1.9 

2.0 

2.2 

2.7 

22 

2 1 

3.4 

4.4 



















23 

1 

13.5 

1 2. 0 

9.3 

7.9 

6.6 

5.5 

4.5 

3.7 

3.0 

2.4 

2.0 

1.5 

1.3 

1.6 

1.4 

1.5 

1.3 

1.5 

1.7 

1.5 

23 

2 1 

1.3 

2.6 

3.0 


















24 

1 

14. 1 

12. 5 

10.5 

7.9 

0.5 

5.3 

4.3 

3.5 

2.8 

2.3 

1.8 

1 .7 

2. 1 

1 « a 

1.3 

1*6 

1.5 

1.4 

1.8 

1 .8 

24 

2 1 

2.2 

2. 1 

2 .6 

3.5 

















25 

1 

14.0 

11.4 

9.7 

7.9 

6.5 

5. 4 

4.4 

3.6 

3*0 

2.4 

1.7 

1*4 

1.8 

1.3 

1.9 

1.3 

1.6 

1.5 

1.4 

1.8 

2 5 

2 1 

1*4 

2. 3 

1 .5 

1.9 

3.2 
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1 1 ) 

1 2 ) 

1 3) 

i 4 > 

1 5 ) 

( 6 > 

( 7 ) 

1 5 ) 

( 9 ) 

( 10 ) 

Ul) 

( 12 ) 

1 13) 

1 14) 

( 15 ) 

( 16 ) 

( 17 ) 

( Id ) 

( 19 ) 

( 20 ) 

z 

1 

15.7 

15 . 7 



















J 

1 

15.7 

15.7 

15.7 


















4 

1 

15.7 

15 . 7 

15.7 

15.7 

















5 

1 

15 . 7 

15 . 7 

15 . 7 

15.7 

15.7 
















6 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 















7 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 














d 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

IS . 7 

15.7 

15.7 













9 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 

15.7 












10 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 

15.7 

15.7 











1 1 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15. 7 

15.7 

15.7 

15 . 7 

15.7 










12 

1 

1 5 . r 

15 # 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 

15*7 









1 J 

1 

15.7 

15 # 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15*7 

15,7 

15.7 








14 

1 

15 . 7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15«7 

15.7 







15 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 






16 

1 

15.7 

15 . 7 

15 .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 





1 7 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 




la 

1 

Id . 7 

15 . 7 

15.7 

15.7 

15 .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 

15.7 

15 . 7 

15.7 

15.7 

15.7 



1 9 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 

15.7 

15 . 7 

15.7 

1 5.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15.7 


20 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

1 5.7 

15.7 

15.7 

15.7 

15.7 

1 5.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

21 

1 

15.7 

1 5 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15*7 

15.7 

16.7 

15.7 

15.7 

15.7 

15.7 

15.7 

2 1 

2 1 

15.7 




















22 

1 

15 . 7 

15 . 7 

15.7 

15.7 

Id. 7 

15 . 7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

22 

21 

15.7 

15 . 7 



















2 J 

1 

15.7 

15 . 7 

15 .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

0.5 

0.5 

0 . 1 

0 . 1 

o 

• 

o 

0*0 

0 . 1 

0.1 

0.5 

23 

21 

0.5 

15.7 

15.7 


















24 

1 

15.7 

15 . 7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

0.2 

0.2 

0.0 

0*0 

0.0 

0.0 

15.7 

24 

21 

15.7 

15 . 7 

15.7 

15.7 

















25 

1 

15.7 

15 . 7 

15 . 7 

15 . 7 

15.7 

15.7 

15 . 7 

1 5.7 

15 . 7 

15.7 

15.7 

15.7 

0 * 1 

0. 1 

0.0 

0*0 

0.0 

0.0 

0.0 

0*0 

25 

21 

0.0 

0. U 

0 .o 

0.6 

15.7 


















Tables - f= .0001 





MATRIX 

ELEMENTS : 

= - DLOG 10 ( lEXACT - COMPUTEDl / EXACT ) ; 

REAL PARTS OF 

POLYRT COMPUTED ROUTS 



( 1 ) 

( 2 > 

{ 3 > 

1 4 ) 

( 5 ) 

( 61 

( 7 ) 

( 6 ) 

( 91 

( 10 ) 

ill ) 

( 12 ) 

( 13 ) 

( 14 ) 

( 15 ) 

( 16 ) 

( 17 ) 

( 18 ) 

i 19 ) 

( 20 ) 

4 

1 

15.9 

15.9 

lb . 9 

15.9 

















6 

1 

1 6 . 2 

16.2 

11.4 

11.4 

11.5 

11*5 















6 

1 

lb .7 

15.7 

1 1 .9 

11.9 

11.6 

1 1.8 

12.3 

12.3 













10 

1 

15.5 

15.5 

12.2 

12.2 

12*0 

12.0 

12.3 

12.3 

12.9 

12.9 











12 

1 

1 b .6 

15.6 

12.4 

12.4 

12.1 

12.1 

12.2 

12.2 

12.7 

12.7 

13.5 

13*5 









14 

1 

14.8 

14.6 

12.5 

12.5 

12.2 

12.2 

12.3 

12.3 

12*5 

12.5 

11.0 

11.0 

11.1 

11.1 







16 

1 

Ib.b 

lb . 5 

12 .6 

12*6 

12.3 

12.3 

12*3 

12.3 

12.6 

12.6 

9.2 

9.2 

9*0 

9.0 

9.4 

9.4 





16 

1 

lb . 1 

15 . 1 

12.7 

12.7 

12.3 

12.3 

12.2 

12.2 

12.5 

12.5 

9.0 

9.0 

8.7 

8.7 

8.8 

8.8 

9.4 

9.4 



20 

1 

15.1 

15*1 

12.7 

12.7 

12.3 

12.3 

12.4 

12.4 

12.9 

12*9 

9.2 

9.2 

8.6 

8.8 

8.9 

8.9 

9.4 

9.4 

10.6 

10.6 

22 

22 

1 

2 1 

15.3 

10*9 

15.3 

10.9 

12.7 

12.7 

12.3 

12.3 

12.1 

12.1 

11.7 

11.7 

9.3 

9.3 

6*9 

8.9 

8.8 

8.6 

9.1 

9 . 1 

9*7 

9.7 

24 

1 

15.1 

15 . 1 

12.6 

12.6 

12.4 

12.4 

12 . 1 

12.1 

11.5 

11*5 

9.5 

9.5 

8.4 

6*4 

6.0 

8.0 

7.9 

7.9 

8.0 

8.0 

24 

2 1 

6.4 

6 . 4 

9.2 

9*2 













26 

1 

15.1 

15 . 1 

12.6 

12.8 

12.4 

12.4 

12.4 

12.4 

12*2 

12.2 

9.6 

9.6 

8.0 

6.0 

7.5 

7.5 

7.2 

7.2 

7.3 

7.3 

26 

2 1 

7.0 

7.6 

6.1 

6.1 

9.0 

9.0 















26 

1 

15 . J 

15.3 

12.6 

12.6 

12.4 

12.4 

12.7 

12.7 

11.2 

11.2 

9.6 

9.6 

7.8 

7.8 

7.1 

7.1 

6.8 

6.8 

6.8 

6.8 

28 

21 

7.0 

7 . 0 

7.7 

7.7 

6.0 

8.0 

6.3 

8.3 








30 

1 

15.2 

15.2 

12.9 

12.9 

12.6 

12.6 

12.2 

12.2 

11.1 

1 1.1 

9.6 

9.6 

7.6 

7.6 

6.9 

6.9 

6*6 

6*6 

6.5 

6.5 

30 

2 1 

6.6 

6 . 6 

6.6 

6.6 

7.3 

7.3 

7.9 

7.9 

8*6 

6.6 








32 

1 

15.4 

15.4 

12 .6 

12.8 

12.4 

12.4 

12.7 

12.7 

11.2 

1 1.2 

9.6 

9.6 

7.5 

7.5 

6.7 

6.7 

6.4 

6*4 

6.2 

6.2 

32 

21 

6.2 

6.2 

6.4 

6.4 

6.6 

6.8 

7.3 

7.3 

6.0 

6.0 

8.9 

8.9 









34 

1 

lb . 2 

15.2 

12.6 

12*6 

12.5 

12.5 

12.1 

12.1 

1 1.0 

11.0 

9.7 

9.7 

7.4 

7.4 

6*6 

6.6 

6.2 

6*2 

6.0 

6.0 

34 

2 1 

6.0 

6.0 

6.1 

6 * 1 

6.6 

6.6 

7.0 

7.0 

6.8 

6 . 8 

7.1 

7.1 

7.9 

7.9 






36 

1 

15.1 

15.1 

12.6 

12.6 

12.3 

12.3 

12.8 

1 2.8 

11.1 

11*1 

9.8 

9,6 

7.2 

7.3 

6.5 

6 . S 

6.1 

6.1 

5*8 

5.8 

36 

2 1 

5.6 

5.8 

5.9 

5.9 

6.1 

6 . 1 

6.5 

6.5 

7.7 

7.7 

7.2 

7.2 

7.4 

7.4 

8.2 

8.2 



38 

1 

Ib .2 

lb .2 

12.9 

12.9 

13.2 

13.2 

11.6 

11.6 

11.3 

11.3 

10.1 

10.1 

7.2 

7.2 

6.4 

6.4 

6.0 

6.0 

5.7 

5.7 

36 

2 1 

5.7 

5.7 

7 .2 

7.2 

5.4 

5.4 

5 . 1 

5 . 1 

5.1 

5 . 1 

5*4 

5.4 

5.8 

5.8 

6.4 

6.4 

7.5 

7.5 




Table 4 - f = .0001 





MATRIX 

tLEMENTS 

= -DL0G10( (EXACT - COMPUTED ) /EXACT ) 

; IMAGINARY PARTS OF POLYRT COMPUTtO ROOT' 



( 1 ) 

( 2 ) 

( 3) 

( 4) 

< 51 

( 6) 

( 7) 

( dJ 

( 9) 

( lOJ 

(11) 

( 12) 

( 13) 

( 14) 

( 15) 

( 16) 

(17) 

(16) 

(19) 

(20) 

4 

1 

15.7 

15. 7 

15.7 

15.7 

















6 

1 

15.7 

15. 7 

15.7 

15.7 

15.7 

15.7 















a 

1 

15.7 

15.7 

16.0 

16.0 

15.5 

15.5 

15.7 

15.7 













10 

1 

15.7 

15.7 

15.7 

15.7 

15.1 

15.1 

15.6 

1 5.6 

15.9 

15.9 











12 

1 

15.7 

15.7 

15.3 

15.3 

14.9 

14.9 

15.7 

15.7 

14.3 

14.8 

15. 1 

15.1 









14 

1 

15.7 

15.7 

15.2 

15.2 

14.5 

14.5 

14.2 

14.2 

14.0 

14.0 

14.0 

14.0 

14.5 

14.5 







16 

1 

15.7 

15. 7 

14.8 

14.6 

14.5 

14.5 

14.5 

14.5 

13.8 

13. d 

13.4 

13.4 

13.4 

13.4 

13.8 

1 3.8 





la 

1 

15.7 

15.7 

14. d 

14.8 

13.9 

13.9 

13.3 

13.3 

12.9 

12.9 

12.6 

12.6 

12.6 

12.6 

12.8 

12.6 

13.5 

13.5 



20 

1 

15.4 

15.4 

15.7 

15.7 

14.0 

14.0 

13.3 

13.3 

13.1 

13«1 

13*4 

13»4 

12.3 

12*3 

12.1 

12.1 

12.3 

12.3 

12.9 

12.9 

22 

1 

15.2 

15.2 

14.7 

14,7 

14.7 

14.7 

13.7 

13.7 

12. a 

12.6 

12.2 

12.2 

12.0 

12.0 

12.0 

12.0 

12.1 

12.1 

12.4 

12.4 

22 

21 

13.0 

13. 0 


















24 

1 

15.2 

15. 2 

14.3 

14. d 

14.2 

14.2 

13.6 

13.6 

12.5 

12.5 

11.9 

11.9 

11.6 

11.6 

1 1.5 

11.5 

11.6 

1 1.6 

1 1 .8 

1 1 .8 

24 

21 

12*2 

12.2 

13.1 

13.1 

















26 

1 

15.4 

15.4 

14.7 

14.7 

14.5 

14.5 

13.1 

13.1 

12.4 

12*4 

1 1 .d 

1 l.d 

1 1.5 

11.5 

1 1.6 

11.6 

1 1.9 

11.9 

11.1 

11.1 

26 

2 1 

10.9 

10.9 

11.2 

11.2 

1 1 .8 

1 l.d 













2d 

1 

15.4 

1 5. 4 

15.0 

15.0 

13.7 

13.7 

12.9 

12.9 

12.2 

12.2 

11.6 

11.6 

11.2 

11.2 

lO.d 

10.3 

10.4 

10*4 

10.2 

10.2 

28 

2 1 

10.2 

10.2 

10.4 

10.4 

10.8 

lO.d 

11.6 

11.6 













JO 

1 

15.4 

15.4 

15.2 

15.2 

13.9 

13. 9 

13.5 

13.5 

13.2 

13.2 

12.5 

12.5 

11*3 

11.3 

10.5 

10.5 

10.0 

10.0 

9.7 

9.7 

30 

2 1 

9.t 

9.6 

9.7 

9. ^ 

10.0 

10.0 

10.4 

1 0.4 

11*2 

11.2 











32 

1 

15.4 

15. 4 

14.6 

1 4 . 6 

14.1 

14.1 

13.6 

13.6 

12.3 

12.3 

11.6 

1 1*6 

11.1 

11.1 

10.7 

10.7 

10.3 

10.3 

9*8 

9.8 

32 

2 1 

9.5 

9. 5 

9.4 

9.4 

9.5 

9.5 

9.9 

9.9 

10.6 

10*6 

12.8 

12. d 









34 

1 

15.4 

15.4 

14.7 

14.7 

14.1 

14.1 

13.4 

13.4 

12.3 

12.3 

11. d 

1 l.d 

11.5 

11.5 

10.3 

10.3 

9.7 

9.7 

9.4 

9.4 

34 

2 1 

9. 3 

9.3 

9.4 

9.4 

9.9 

9.9 

10.0 

10.0 

9.9 

9. 9 

l(U2 

10.2 

11.0 

11.0 







36 

1 

15.7 

15.7 

1 4.5 

14.5 

13.6 

13.6 

13.9 

13.9 

12.4 

12.4 

14.5 

14.5 

1 O.d 

10.8 

9.9 

9.9 

9.4 

9.4 

9.0 

9.0 

36 

21 

d. d 

d. 3 

d .d 

d .6 

9.0 

9.0 

9.3 

9.3 

9.9 

9.9 

10.6 

10.6 

11.2 

1 1.2 

12.0 

12.0 





3d 

1 

15.7 

15.7 

14.4 

14.4 

13.5 

13.5 

13.0 

13.0 

12.8 

12.0 

12.2 

12*2 

I 1 . 1 

11.1 

10.6 

1 0.6 

10.9 

10.9 

9.4 

9.4 

3d 

2 1 

d. ti 

d. 6 

d .6 

8.6 

8.5 

8.5 

d. d 

8.8 

9.5 

9. 5 

a. 9 

8.9 

9.0 

9.0 

9.4 

9.4 

10.3 

10.3 





Table 5 - ? = .9 





MATW IX 

Elements 

= -DLOG10( (EXACT 

- CUMPUTED)/EXACT) ; 

kCAL FAKTS OF 

POLYRT COMPUTED ROOTS 



( 1) 

( 2 ) 

1 3) 

( 4) 

( 5) 

( 6) 

( 7) 

( 8) 

( 9) 

(10) 

(11) 

( 12) 

( 13) 

(14) 

( 15) 

(16) 

( 17) 

(18) 

(19) 

(20) 

4 

1 

15.9 

15. 9 

15.9 

15.9 

















6 

1 

15.1 

15. 1 

14.7 

14.7 

14.8 

14.8 















6 

1 

15.3 

15. 3 

13.6 

13.8 

13.5 

13.5 

13.4 

13.4 













10 

1 

14.7 

14. 7 

13.3 

13.3 

12.3 

12.2 

12.4 

12.4 

13.5 

13.5 











12 

1 

14.2 

14.2 

12.9 

12.9 

1 1 .6 

11.6 

11.5 

1 l.b 

11.3 

1 1.3 

12.1 

12.1 









14 

1 

14.1 

14.1 

12.7 

12.7 

1 1.2 

11.2 

11.3 

1 1.3 

11.0 

11.0 

10.9 

10.9 

11.7 

11.7 







16 

1 

1 

13.9 

1 1 .6 

11.6 

10.6 

10.6 

9.5 

9.5 

9.2 

9.2 

9. 1 

9.1 

10.1 

10.1 

10.4 

10.4 





la 

1 

13.2 

13.2 

1 1 .7 

11.7 

9 .7 

9.7 

9.5 

9.5 

8.4 

8.4 

9.0 

9.0 

8.5 

8.5 

9.5 

9.5 

9.6 

9.6 



20 

1 

13.7 

13.7 

11.9 

11.9 

9.9 

9.9 

6.5 

8.5 

8.2 

8.2 

7.2 

7.2 

7.6 

7.6 

7.4 

7.4 

7.7 

7.7 

8.8 

8.6 

22 

1 

13.5 

13. 6 

10.5 

10.5 

9.2 

9.2 

a.4 

8.4 

7.6 

7.6 

6.7 

6.7 

7.2 

7.2 

6.4 

6.4 

6.8 

6.8 

7.3 

7.8 

22 

21 

8.7 

a. 7 



















24 

1 

13.4 

13.4 

10.4 

10.4 

8.7 

6.7 

7.4 

7.4 

8.6 

6.6 

6. 1 

6.1 

5.6 

5.6 

6. 1 

6. 1 

5.6 

5.6 

6.0 

6.0 

24 

2 1 

8.7 

6. 7 

7.1 

7.1 

















26 

1 

12.9 

12.9 

10.0 

10.0 

8.3 

a. 3 

8. a 

6.8 

5.8 

5.6 

5.2 

5.2 

4.7 

4.7 

4.2 

4.2 

4.6 

4.6 

4.7 

4.7 

26 

21 

4.7 

4. 7 

5 .3 

5.3 

6.5 

6.5 















28 

1 

13. 2 

13. 2 

10.0 

10.0 

8.2 

6.2 

6. 8 

6.8 

5.5 

5.5 

5.3 

5. 3 

4.2 

4.2 

4.5 

4.5 

4. 1 

4.1 

3.9 

3.9 

28 

21 

4.1 

4. 1 

4 .5 

4.5 

5.4 

5.4 

7.0 

7.0 













30 

1 

12.8 

12 . a 

9.8 

9.6 

7 .9 

7.9 

6.3 

6.3 

5.1 

5. 1 

4.1 

4.1 

3.9 

3.9 

2.8 

2.8 

2.7 

2.7 

2.5 

2.5 

30 

2 1 

2.3 

2. 3 

2./ 

2.7 

3.5 

3.5 

J . 6 

3.6 

4.5 

4.5 











32 

1 

12.8 

12. o 

9.5 

9. 5 

7.6 

7.6 

7.9 

7.9 

4.6 

4.6 

3.7 

3.7 

2.9 

2.9 

2.4 

2.4 

2.0 

2.0 

2.2 

2.2 

32 

21 

1.7 

1 . 7 

2.3 

2.3 

2.4 

2.4 

2.8 

2.8 

4.1 

4.1 

4.8 

4.8 









34 

1 

14.0 

1 4. 0 

9 .5 

9.5 

7.3 

7.3 

5.6 

5.6 

8.0 

6.0 

3.2 

3.2 

3.0 

3.0 

2.3 

2.3 

1.4 

1.4 

1.7 

1.7 

34 

2 1 

1.3 

1.3 

1 .2 

1.2 

2.3 

2.3 

1.1 

1.1 

1.7 

1.7 

2.1 

2.1 

3.7 

3.7 







36 

1 

12«b 

12.8 

10.7 

10.7 

7.1 

7. 1 

5.5 

5.5 

4.2 

4.2 

3. 1 

3.1 

3. 1 

3.1 

1.9 

1.9 

1.2 

1.2 

1.3 

1.3 

36 

21 

0.9 

0.9 

1.0 

1.0 

1 .4 

1.4 

1.8 

l.b 

1.2 

1.2 

2.6 

2.6 

1.8 

1.8 

2.6 

2.6 







Table 6 - ? = .9 





MATKIX 

t £ Lt ,Mt,NT 3 

= - 0 L 0 G 10 ( <EXAC 



( 1 ) 

( 2 ) 

i 3 ) 

( 4 ) 

( 5 ) 

C 6 ) 

( 7 ) 

( 8 ) 

4 

1 

lb «3 

lb . 3 

15.4 

15.4 





6 

1 

14.7 

14.7 

14.3 

14.3 

14.9 

14.9 



a 

1 

1 * 4.4 

14.4 

13.4 

13.4 

12.9 

12.9 

14 . 1 

14 . 1 

10 

1 

15.5 

15.5 

12.6 

12.6 

12.4 

12.4 

12.0 

12.0 

12 

1 

14.5 

14.5 

12.2 

12.2 

11.5 

1 1.5 

10.9 

10.9 

14 

1 

13.3 

13.3 

1 1 .8 

1 1.8 

12.0 

12.0 

10.7 

10.7 

16 

1 

13.4 

13.4 

1 1 .4 

11.4 

9.9 

9.9 

9.2 

9.2 

la 

1 

13.6 

13.6 

10.9 

10.9 

9.8 

9.8 

8.5 

6.5 

20 

1 

13.6 

13.6 

1 1 .2 

11.2 

9.5 

9.5 

8.5 

8.5 

22 

1 

13.0 

13.0 

11 .8 

1 l.d 

9.1 

9.1 

8.4 

8.4 

22 

21 

a.o 

a . 0 







24 

1 

13.2 

13.2 

10.4 

10.4 

8 . a 

e.e 

9.0 

9.0 

24 

21 

b .9 

5.9 

7.1 

7.1 





26 

1 

12.4 

12.4 

9 .9 

9.9 

8.1 

e . 1 

6.8 

6.3 

26 

21 

4.6 

4.8 

b .b 

5.5 

6.1 

6 . 1 



28 

1 

12.2 

12.2 

9.6 

9.6 

7.9 

7.9 

c . 5 

6.5 

2 d 

21 

4.6 

4.6 

4.3 

4.3 

4.7 

4.7 

5.6 

5 . 6 

30 

1 

12.2 

12.2 

10.1 

10.1 

7.4 

7.4 

5. .9 

5.9 

30 

21 

3.2 

3.2 

2 .2 

2.2 

2.5 

2.5 

3.4 

3.4 

32 

1 

12.0 

12 . 0 

10.0 

10.0 

7.3 

7.3 

5.5 

5.5 

32 

21 

1 .a 

1 . 8 

1 .5 

1.5 

2.0 

2.0 

3.0 

3.0 

34 

1 

12.0 

12 . 0 

9.6 

9.6 

7.9 

7.9 

5 . 4 

5.4 

34 

2 1 

0.3 

0.3 

0.9 

0.9 

0.8 

0.8 

0.9 

0.9 

36 

1 

12.1 

12 . 1 

10.4 

10.4 

7.2 

7.2 

4.9 

4.9 

36 

2 1 

0.3 

0.3 

0.9 

0.9 

0.5 

0.8 

0.7 

0.7 


- COMPUTeO)/EXACT) ; IMAGINARY PARTS OF POLVRT COMPUTED ROOT 


( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 

( 13 ) 

( 14 ) 

( 15 ) 

( 16 > 

( 17 ) 

( 18 ) 

( 19 ) 

( 20 ) 

12.4 

12.4 











12.2 

12.2 

11.7 

11.7 









10.7 

10.7 

11.3 

11.3 

1 1.3 

11.3 







6.8 

8.8 

9.2 

9.2 

9 . 1 

9 . 1 

9.9 

9.9 





8.7 

8.7 

8.0 

a.o 

9.2 

9.2 

8.5 

8.5 

9.2 

9.2 



7.4 

7.4 

8.0 

8.0 

6.8 

6.8 

7.0 

7.0 

9.6 

9.6 

8.4 

8.4 

7.0 

7.0 

6.7 

6.7 

6.1 

6.1 

6.6 

6.6 

6.7 

6.7 

7.0 

7.0 

6.3 

6.3 

5.6 

5.6 

5.7 

5.7 

5.2 

5.2 

5.8 

5.8 

5.5 

5.5 

5.6 

5.6 

4.7 

4.7 

4.3 

4.3 

4.5 

4.5 

3.9 

3.9 

4.0 

4.0 

7.6 

7.6 

4.4 

4.4 

4.6 

4.6 

3.6 

3.6 

3.6 

3.6 

4.1 

4. 1 

4.8 

b.O 

4.8 

5.0 

4 . 0 

4.0 

3.0 

3.0 

3.4 

3.4 

2.2 

2.2 

2.0 

2.0 

4.5 

3.3 

4.5 
3 . 3 

3.4 

4.6 

3.4 

4.6 

2.6 

2.6 

2.0 

2.0 

2.0 

2.0 

1.5 

1.5 

3.3 

1.2 

3.0 

1 . 2 

3.5 

1.9 

3.5 

1 .9 

2.3 

2.6 

2.3 

2.6 

1.5 

1.5 

1.4 

1.4 

0.8 

0.8 

3.6 

0.7 

3.6 

0.7 

3.5 

1.1 

3.5 

1.1 

4 . e 
2.2 

4.6 

2.2 

1.2 

3.3 

1.2 

3.3 

2 . 1 

2.1 

0.9 

0.9 



Table 7 - f = 1.0 


MATKIX aLcMt-iNTS = -DLUolO( (dXACr - COMPOTtO)/EXACT) ; HEAL PAWTS OF POLVRT COMPUTED ROuTS 




( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 

( 7) 

( 8 ) 

( 9 ) 

( 10 ) 

(ID 

( 12 > < 

( 13 ) 

( 14 ) 

< 15 ) 

C16 ) 

1 17 ) 

( 18 ) 

( 19 ) 

( 20 ) 

4 

1 

8.3 

8.3 

lb ./ 

15 . / 

















6 

1 

8.3 

8 . o 

7.9 

7.9 

7.3 

/.3 















8 

1 

8.4 

8 . 4 

7.-3 

7.3 

7.7 

7.7 

7.8 

7.8 













10 

1 

8.4 

8 . 4 

6.2 

6.2 

b.5 

5.5 

10.3 

1 0.3 

10.5 

10.5 











12 

1 

8.4 

8.4 

b • 4 

5.4 

5.1 

b . 1 

5.1 

b . 1 

4.8 

4.6 

5.0 

5.0 









14 

1 

8.5 

8.5 

5*4 

b.4 

8.0 

8.6 

7.7 

7.7 

8.1 

8. 1 

7.8 

7.8 

8.b 

8.5 







1 b 

1 

8»b 

8 . D 

b .4 

5.4 

4.3 

4.3 

3.7 

3.7 

3.3 

3.3 

3.3 

3.3 

3.5 

3.5 

4 . 1 

4 . 1 





18 

1 

8.b 

8 . 5 

b .4 

b.4 

4 .2 

4.2 

3.4 

3.4 

3.0 

3.0 

2.9 

2.9 

5.2 

5.2 

5.8 

5.8 

6 . 1 

6.1 



20 

1 

8.6 

8.5 

5.4 

b.4 

3.7 

3 . 7 

2.9 

2.9 

2.5 

2.4 

2.5 

2*5 

3.6 

3.6 

4.0 

4.0 

3.8 

3.6 

4.9 

4.9 

22 

1 

8.5 

8 . 5 

5 .4 

S.4 

3.8 

3.8 

2.8 

2.8 

2.3 

2.2 

3.0 

3.0 

2.4 

2.4 

2.8 

2.8 

2.5 

2.5 

2.7 

2.7 

22 

2 1 

3.6 

3.8 



















24 

1 

8.5 

8 . 5 

5.3 

5.3 

3.1 

3 . 1 

2.2 

2.1 

1.5 

1.1 

1.0 

1.3 

1.0 

1.3 

1.1 

1.3 

1.2 

1.2 

1.4 

1.1 

24 

21 

1.6 

1.2 

1.7 

1.9 

















26 

1 

8.6 

8.6 

5.3 

b.3 

3.0 

3.0 

2.1 

2.0 

1.5 

1.2 

1.4 

1.4 

1.2 

1.2 

1.4 

1 .4 

2.0 

2.0 

1.7 

1.7 

26 

2 1 

1.4 

1.4 

1 .5 

1.5 

2.1 

2 . 1 















28 

1 

8.6 

8 . 6 

b .4 

5.4 

2 .8 

2.7 

1.7 

1.4 

1.0 

1.7 

0.8 

2.1 

0.8 

2.5 

0.9 

1.7 

1.0 

1.3 

1.3 

1.1 

28 

21 

1.7 

1.0 

2.4 

1.0 

1.9 

1.1 

2*8 

1.8 













30 

1 

8.6 

8.6 

5.2 

5.2 

2.5 

2.5 

1.5 

1.2 

0.6 

2.3 

0.8 

1.5 

0.8 

1.6 

o.a 

4.2 

1.0 

1.2 

1.2 

1.2 

30 

21 

1 .V 

1 . ^ 

1 .4 

1.4 

1 .1 

1 . 1 

1*1 

1.1 

1.3 

1.3 











J2 

1 

8.0 

8.3 

b .2 

b.2 

2.5 

2.5 

1.4 

1.3 

0.8 

1.6 

0 . 7 

1.2 

0.7 

1 . 1 

0.7 

1.2 

0.8 

1.5 

0.9 

3.5 

32 

21 

1.0 

1*9 

1.1 

1.2 

1 .7 

0.9 

1.4 

0.8 

1.1 

0.3 

1 • 1 

1.0 









J4 

1 

8.2 

8.2 

5.1 

5 . 1 

3.7 

3.8 

0.9 

1.3 

0.6 

0.6 

0.6 

0.6 

0 . 7 

0.7 

O.d 

0.8 

0.6 

1.2 

0.8 

2.2 

34 

21 

1.1 

1 . 1 

2*0 

0.7 

1 .2 

1.2 

1.4 

1.4 

1.0 

1.0 

0.9 

0.9 

0.9 

0.9 







36 

1 

8.1 

8 . 1 

5.2 

5.3 

3.2 

3.2 

0.6 

0.7 

O.b 

0.6 

0.4 

0.7 

0.5 

0.8 

0.6 

1.0 

0.7 

1.4 

0.9 

1.7 

36 

21 

1 .u 

1.4 

0.9 

0.9 

1.3 

1.3 

1.5 

1.5 

1.0 

1.0 

0.8 

0.8 

0.7 

0.7 

0.7 

0.7 







Table 8 -f= 1.0 


MAfr<lA £L = Mt;Nlo - -DLOG 10 ( ( EXAC T - COMPUTtiO) /EXACT ) ; IMAuXNAKV PAkTS OP POLYHT CUMPUTrO HOOTS 




( 1 ) 

( 2 ) 

( 3 ) 

( 

( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 

( 0 ) 

( 10 ) 

( 11 ) 

( 12 ) 

( 13 ) 

( 14 ) 

( 15 ) 

( 16 ) 

( 17 ) 

( 16 ) 

( 19 ) 

( 20 ) 

4 

1 

lb . 7 

15.7 

7.4 

7.4 

















6 

1 

lb . 7 

15.7 

15.7 

15.7 

Id . 7 

15 . 7 















8 

1 

lb .7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

3.5 

3.5 













1 U 

1 

16.7 

lb .7 

15.7 

15.7 

15.7 

15.7 

^•9 

4.9 

4^9 

4.9 











12 

1 

Id . 7 

15 . 7 

Id .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 









X 4 

1 

16.7 

15.7 

15.7 

15.7 

4.2 

4.2 

3.0 

3.6 

3.3 

3.3 

3.4 

3.4 

4.0 

4.0 







16 

1 

lb .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15 * 7 

lb .7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 





Id 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

15.7 

3.0 

3.0 

2.2 

2*2 

2.7 

2.7 



20 

1 

15.7 

15.7 

15.7 

15.7 

15.7 

15 . 7 

15 . 7 

1 5 . 7 

15.7 

15.7 

15.7 

15.7 

1 .3 

1.3 

1.2 

1.2 

1 .4 

1 .4 

2.1 

2.1 

22 

22 

1 

2 1 

15.7 

1.6 

15 . 7 
1.6 

15.7 

15.7 

15.7 

15.7 

lb .7 

16.7 

15.7 

15.7 

1*4 

1*4 

0.7 

0.7 

0.5 

0.5 

0.6 

0.6 

0.9 

0.9 

24 

24 

1 

21 

15.7 
0 . 1 

15 . 7 
0.5 

15.7 
0 .5 

15.7 

15.7 

15.7 

15 . / 

15.7 

15.7 

15.7 

1.1 

1.1 

0.2 

0.2 

0.1 

0 . 1 

0 . 0 

0.0 

0.0 

0.0 

0*1 

26 

26 

1 

21 

16.7 
0 . u 

15.7 
0. 0 

15.7 
0 .3 

15.7 

0.3 

15.7 

15.7 

15.7 

15.7 

15 . 7 

1 5 . 7 

15.7 

15.7 

O.d 

0*3 

0 * 1 

0.1 

0.0 

0*0 

0.0 

0.0 

0.0 

0.0 

26 

26 

1 

21 

15.7 

0.0 

15 . 7 
o . 0 

15.7 

0.0 

15.7 

0.0 

15.7 

0.0 

15 . 7 
0.3 

15.7 
0 . 3 

15.7 
1 5 . 7 

15.7 

0.2 

0.2 

0.0 

0.0 

0.0 

0.0 

0.0 

0*0 

0.0 

0.0 

0*0 

60 

30 

1 

21 

15.7 

0.0 

15 . 7 
0 . J 

15.7 

0.0 

15 . 7 
0.0 

15.7 

0.0 

15.7 

0.0 

15.7 

0.0 

0.6 

0.0 

0.0 

0.3 

0 . 1 
0.3 

0 . 1 

0.0 

0*0 

0.0 

0.0 

0*0 

0.0 

15.7 

0.0 

0.0 

32 

32 

1 

21 

15.7 

0.0 

15.7 
0 . 0 

15 .7 
0.0 

15.7 

0.0 

15.7 

0.0 

15.7 

0.0 

15 . 7 
0.0 

0.6 

0.0 

0.6 

0.0 

0.1 

0.0 

0 . 1 
0*0 

0.0 

15.7 

0*0 

0.0 

0.0 

0.0 

0.0 

0.0 

0*0 

0.0 

34 

34 

1 

2 1 

15.7 

0.0 

15 . 7 
0.0 

15.7 

0.0 

15.7 

15.7 

15.7 

0.0 

15.7 

0.0 

15.7 

0.0 

13.7 

0*0 

0.0 

0.0 

0 . 0 
0 . 0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

15.7 

0.0 

0.0 

0.0 

0.0 

0.0 

36 

36 

1 

2 1 

15.7 

0.0 

15.7 
15 . 7 

15.7 
0 .0 

15.7 

0.0 

1.7 

0.0 

1.7 

0.0 

15.7 

0.0 

0.0 
0 . 0 

0*0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0*0 

0.0 

0.0 

0*0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 



Tables 7 and 8 present results associated with equation 6 for f = 1 .0. These results are intended to 
provide data at the limit of the breakdown of computation, that is, when roots are exactly coinci- 
dent in the complex plane, or when damping is critical. 

Data for over-damped systems is not presented, since all roots would then become real, and data 
presented in Tables 1 and 2 would then become relevant to this situation. 

In summary, the data presented shows that polynomials in power series form are adequate for 
low- to moderate-order oscillatory systems which do not have any near-critically damped system 
frequencies. For oscillatory systems with nearly critically damped system frequencies and for 
non-oscillatory systems, a yellow caution flag is raised. The disturbing aspect of the data derived is 
that numerical accuracy degrades gradually. Recognizing that users always attempt to push meth- 
ods a little too far, it is inevitable that if a change in approach is not made, false conclusions will 
inevitably result. 

SYSTEM EQUATIONS, PRIMARY DATA 

Specification of a mathematical simulation model begins at the point at which all input param- 
eters are physically measurable items of data (e.g., length, mass, feedback gain, position, orien- 
tation). The computing problems which lie along the path from this starting point to the need to 
evaluate transfer functions in the complex frequency domain must now be examined. 

The most extensively used representation for a feedback control system is the block diagram. 
The parameters defined in each block of the block diagram, before reduction to canonical form, can 
normally be considered as “primary data,” or physically measurable items. From the block diagram 
it is possible to set up both time and frequency domain equations which characterize system response. 

If the system is continuous, a system of first order ordinary differential equations of the form: 


X = AX + BU 

(8) 

Y = CX 

(9) 


can be directly set up, where X are state variables, Y and U are output and input variables, and the 
coefficient matrices A, B, and C have all elements defined in terms of primary data. If the sys- 
tem is discrete, a system of first order finite difference equations of the form: 


X(k+1) = AX (k) -I- BU (k) 


( 10 ) 
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Y (k) = CX (k) 


( 11 ) 


can be directly set up with all matrix elements defined in terms of primary data. 

Sampled data systems do not lead as readily to continuous or discrete time equations. They can, 
however, be defined in terms of linear mixed operation equations (see Reference 10) with all 
matrix elements defined in terms of primary data. From this set of mixed operation equations it 
is possible to set up, via a straightforward procedure, a system of first order finite difference equa- 
tions of the form given by equations 10 and 11. This set-up procedure, however, requires the 
computation of both a matrix exponential and an integral involving the matrix exponential. A 
detailed study of the condition and of the relevant associated computational operations has been 
done by van Loan, (see References 22, 30, 31, and 32). With respect to computational error, if 
subroutine FADE (Reference 32) is used, then the user can define on input an acceptable bound for 
the perturbation matrix E. The algorithm uses this data to set control parameters which guarantee 
that results are exact for some matrix A + E. With respect to condition, van Loan (see Reference 
30) concludes that “the matrix exponential problem is relatively well-conditioned when A [in e^*] 
is a normal matrix and more poorly-conditioned where A has a defective eigensystem.” These 
results do not answer the important question of whether poorly-conditioned eigensystems imply 
sensitivity of the exponential. This is an area for future research. 

One way to get a “feel” for whether or not A is poorly conditioned with respect to exponentiation 
is to attempt to find the nonsingular transformation matrix <p which will reduce it to block diagonal 
form (see method 18 in Reference 22). This can be done via subroutine BDIAG, developed and 
given in Reference 1 , or by its extended version, namely BLOCK IT which is discussed in Reference 
9. If the resultant blocking is poor, that is if the algorithm yields blocks much greater than one or 
two and <p has a large condition number, with respect to matrix inversion, “beware.” 

Once the feedback control equations are defined in the form given either by equations 8 and 9 or 
equations 10 and 11, numerous methods exist for obtaining transfer functions as polynomial 
ratios in factored polynomial form. In a relative sense, poles are easy to obtain. They are simply 
the eigenvalues of A which are obtainable via application of the EISPACK algorithms (Reference 27). 
Computation of zeros is considerably more difficult, as proven by the profusion of “better” methods 
in the literature, (see References 3, 4, 5, 6, 17, 18, 20, 23 and 26). 

If the objective is the computation of transmission zeros; i.e., in a rough sense those complex 
frequencies at wliich transmission through the system is blocked, then application of the QZ-method 
is currently the best available (see Reference 18). For a definition of the QZ-method, see References 
21 and 29, and for its subroutine code, see Reference 1 1. 

If the objective is the computation of the zeros for particular transfer functions in the transfer 
function matrix, any of the referenced methods can be used if the system is of low order. Caution 
must be exercised if the system is of moderate to large order. Davison, in Reference 5, page 481, 
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cautions against the use of Brockett’s method (Reference 4). Davison’s method (Reference 5) is 
unsatisfactory because it is based upon introducing an arbitrarily large number into a matrix, the 
eigenvalues of which are the sought-after zeros. This approach has the net effect of increasing the 
norm of the matrix and hence the error associated with follow-on eigenvalue computation. Eigen- 
values, computed via the QR method are exact for some matrix A + E where the norm of the 
perturbation matrix E is proportional to the norm of A (see Reference 36, pages 353 and 367). 

Sandberg and So’s method (Reference 26) relies on a correct determination of numerator degree. 
This is a difficult numerical task for moderate- to large-order problems; an incorrect determination 
leads to incorrect results. Bollinger’s method (Reference 3), Marshall’s method (Reference 20), and 
Patel’s method (Reference 23) each rely upon the computation of the roots of a polynomial in 
power series form. This, for moderate- to large-order systems, is an ill-conditioned computing 
problem. Kaufman’s method (Reference 1 7) is based upon solving the generalized eigenvalue 
problem and can be readily implemented by use of the QZ-algorithm rather than Kaufman’s 
algorithm. If the QZ-algorithm is used results are exact for the problem (A + E) X = X (B + F) X 
where bounds in both E and F can be defined. In application, the only problem the author has 
encountered with this method is that eigenvalues known a priori to be identically equal to zero are 
computed as “small non-zero numbers.” For the problem of simulating moderate- to large-order 
low frequency systems these small numbers tend to get mixed up with true small magnitude eigen- 
values. This is a problem with the QZ-algorithm, however, and not with Kaufman’s method. 
Perhaps this problem can be corrected by prefacing the entrance into the QZ-method with a pre- 
balancing step analogous to the prebalancing step associated with the QR-method and implemented 
via subroutine BALANC (see Reference 27). The goal would be to isolate eigenvalues which can be 
determined via zero error row-column permutations. In application, this is extremely useful since 
the eigenvalues which are isolated have zero computational error. Frequently these are the zero fre- 
quencies which, if computed, become troublesome small magnitude complex numbers. This is an 
area of future study. 

It is the author’s opinion that for moderate to large order systems the zero-finding algorithm 
NUMS contained in the DISCOS program (Reference 2, page 65, Vol. I and page B-14I, Vol. II) is 
the best available. It does not have any of the disadvantages cited above and has been successfully 
applied on some very difficult poorly-conditioned problems for which the other methods above 
have been less than satisfactory. An upgraded version of NUMS is currently in the DISCOS program 
being distributed by COSMIC. An equivalent version can be easily created by simply using EISPACK 
(Reference 27) and UNPACK (Reference 8) based subroutines where required. These were not 
available when NUMS was originally coded. 

The preceding remarks pertain when the so-called “state variable” approach is utilized. In general 
application it is often desirable to obtain feedback control transfer functions via block diagram 
algebra (see Reference 7). If this approach is applied, the net effect is to arrive at transfer functions 
defined as ratios of polynomials defined in power series form. These are normally obtained via a 
sequence of algebraic operations with functions of polynomials. As previously stated, finding the 
zeros of a moderate- to large-degree polynomial in power series form is frequently an ill-conditioned 
computing problem and should be avoided. 
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If all blocks of the block diagram can be defined as low-order polynomial ratios in either the Laplace 
variable s or the Z-transform variable z then it is possible to write directly from the block diagram a 
higher order state variable equation of the form: 


P(s)Y(s) = BU(s) 


( 12 ) 


where the elements of the square matrix P (s) are ratios of low-degree polynomials in power series 
form and U (s) and Y (s) are input and output variables. 

If the feedback system is sampled data, then a mixed operation higher order state variable equation 
can be written of the form: 


P,l(s) Pj 2 (z) 

1 

> 

1 Y, W 



(S) P^2 (Z) 

> 

1 

[y, (z,I 

1' 

1 

o 



(13) 


where again matrix elements are ratios of low-degree polynomials. In this case input must be 
assumed constant over the sampling interval, hence U (z) and not U (s). 

At present, there exists no general-purpose algorithm to go directly from equation 13 to an equation 
of the form: 


P(z)Y(z) = BU(z) 


(14) 


However, it appears that a combination of the ideas presented in References 10 and 38 will lead to 
the desired transformation. This is an area of future research. 

In order to obtain the elements of the transfer function matrix G (z) or G (s) an inverse of a matrix 
of polynomials in s or z must be computed. That is: 


G (z) = P (z)-» 


(15) 
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On tlic bottom line this calls for the need to compute the zeros of a series of determinantal equa- 
tions, in which matrix elements are low-order polynomials in power series form. It should be noted 
that, at this point, polynomial coefficients can still be considered as “primary data.” 

Let the matrix of rational polynomials, for which the determinant must be computed, be given by: 


H(z) 
d (z) 


(16) 


where d (z) is the polynomial which is the least common denominator and H (z) is a square matrix 
of low-degree polynomials in power series form. If r is the degree of the highest degree polynomial 
in H (z), then it is rather straightforward to write: 


det [H (z)] = det [z"^ + z^‘‘ j +**•+€„] = 0 


(17) 


where the order of each of the matrices Cj , • ♦ •, q equals N, the order of H (z). Moler and 
Stewart (Reference 21) point out that this is the “X-matrix” problem. On page 242 they show that 
it is solvable via the QZ-algorithm. Their approach, however, demands that a generalized eigen- 
value problem of order N(r + 1) be solved. In application, it is frequently possible to set up a 
generalized eigenvalue problem of considerably lower order, as follows. 

Step 1. Scan each column of the matrix H (z) and record the degree of the highest degree poly- 
nomial. Let: 


r. - degree of liighest degree polynomial in column i; 
if r. = 0, reset it equal to 1 . 


Step 2. View z as a differential operator and consider the set of transform equations: 


Hj, (z) - • -H,^(z) 


Hn,(z)---H^n(^)J 


X, (z) 


.X^ (z) 



(18) 
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where each of the state variables in X (z) are higher order state variables. This equation can be 
rewritten as a set of first order state variable equation of the form: 


A X (z) + zBX (z) = 0 


(19) 


where the order of system of equations is R = r^ + r^ + • * * + rj^ . This algorithm is best illustrated 
by a simple example. Consider the system: 


Let: 


"a 3 z 3 + z 2 + 0 b^ z^ + b„ 


'x, (z)' 

1 


C 4 z"* + C 3 z^ do z 


X, (z) j 

= < 

A 

.ejZ^+CiZ + eo /„ z^ + 

1 

,X 3 (Z)J 


[0) 


( 20 ) 



(z) = 

X, 

(z) 

^3 1 

(z) = 

X3 (z) 

^12 

(z) = 

zX 

11 (z) 

^3 2 

(z) = 


(z) 


(z) = 

zX 

12 (Z) 

X33 

(z) = 

ZX33 

(z) 

^14 

(z) = 

zX 

13 (Z) 

X34 

(z) = 

ZX33 

(z) 





X35 

(z) = 

^^3 4 

(Z) 

^2 1 

(z) = 

^2 

(z) 






( 21 ) 
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Directly substitute into equation 32 to obtain: 



( 22 ) 


The important point to note here is that we have arrived at the need for solving a generalized 
eigenvalue problem via a trivial zero-computational-error operation. 

The eigenvalues of equation 22, which is of the form of equation 19, are the values of z which 
satisfy the determinantal equation 17: 


det [H (z)l = 0 


(17) 


These can then be used to rewrite this equation in the desired factor polynomial form: 


det [H(z)l = G*(z-z,)(z-Z 2 ) • • • (z-z^) = 0 


(23) 
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Where M < R and where G* is obtained by evaluating both sides of equation 23 for any value of z 
not equal to an eigenvalue and solving directly for G*. 

M is determined by examining the computed values of z.; j = 1, 2, • • *, R. M is equal to the num- 
ber of computed eigenvalues “not” at numerical infinity in the complex plane (see Reference 21, 
page 241). 

REDUCED ORDER SYSTEM EQUATIONS 

In References 1 and 9, it is shown that a computational algorithm has been developed which can be 
used to compute the transformation matrix (j) which will reduce a real nonsymmetric matrix A to 
block diagonal form. This can be expressed as: 


0-* A(p = diag(Bj, Bj, • • •, B^) 


(24) 


where k is the number of blocks. The transformation matrix is guaranteed to be nonsingular, and 
all elements of it are real. Since each block is quasi-upper triangular, usually 1 X 1 or 2 X 2, each 
block and their associated columns of <p can be interpreted as system frequencies and quasi-system 
modes. In brief consider the problem: 


X = AX + BU 


(25) 


Find <f> and define the coordinate transformation: 


X = 0q 


(26) 


Direct substitution leads to: 


q = A0q + 0-‘ BU 


(27) 
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or by equation 24. 


q = diag (B^ , , • • • , B^ ) q + BU 


(28) 


Transform the equations and obtain: 


q(s) = [s{]-diag(Bj, • • •,Bj^)] ‘ 0 ‘ BU (s) 


(29) 


where LI is the unit diagonal of appropriate dimension. 
Then back transform and obtain 


X (s) = 0 diag [(s)]- Bj)'S (sLl- B 2 )■^ • • *, (st]- B^)'* ] BU (s) (30) 


where partial inversion of the block diagonal matrix has been carried out. 

Making use of the fact that the elements of 0 define system mode observability and that the ele- 
ments of 0‘* define system mode controllability, it is possible to utilize engineering judgement to 
write reduced order transfer functions as partial fraction expansions of the form 


G(s) = 


k N.(s) 


(31) 


where L < K. 

Normally, the degree of all D^. (s) is one or two with the degree of (s) being less than or equal to the 
degree of (s). Larger order blocks lead to larger degree polynomials, which can be determined via 
the “method of Leverrier” presented first in 1840 (Reference 14, page 166). The method is outlined in 
Appendix A. When using the method of Leverrier it is essential to bear in mind that this method is 
not recommended for “large matrices,” since it requires too many operations and tends to be 
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numerically unstable (Reference 14 , page 149 ). The method of Leverrier is applicable, however, for 
the determination of the matrix of polynomials defined by: 


(sd-Bp-‘ 


due to the fact that each matrix B. is of low order. If the matrix B. is moderate to large order, 
blocking is by definition “poor.” This normally implies that the basic computing problem is ill- 
conditioned. This does not imply that the physical system is ill-conditioned. It does, however, imply 
that the analyst, through poor engineering judgment, has set up a set of system dynamics equations 
which defy solution. A few common examples of poor engineering judgment are: 

• An attempt to model system response frequencies which differ by many orders of magnitude; 

• A poor choice of physical units. Ideally, all numbers of significance in the computing prob- 
lem should be of comparable order; and 

• A poor choice of state variables such as, for example, using a Euler angle sequence which is 
subject to “gimbal lock,” using relative rather than absolute coordinates for rate and/or 
displacement measures or vice versa. 

When blocking is poor, “stop.” The analyst is advised to step back and to reexamine the basic 
formulation which led to the troublesome computing problem. It is usually a much simpler prob- 
lem to reduce frequency spectrum, change units, or use a different set of state variables than to 
fight an ill-conditioned computing problem. 

Again, the objective is to write the transfer function G (s) given in equation 31 as a truncated 
partial fraction expansion as a ratio of polynomials both given in factored polynomial form. The 
denominator is the least common denominator; since all B. in equation 29 are quasi-upper triangular 
it is simply the root product of all system frequencies not truncated out. The computation of the 
numerator polynomial, in factored polynomial form, can be computed via the QZ-algorithm by 
recognition of the fact that: 


0 (s) D3 (s) 

Dj (s) 0 D3 (s) 

D, (s) (s) 0 


(s) 1 

Dl (s) 1 

Dl (s) 1 


( 32 ) 


Dj (s) D3 (s) D3 (s) 0 1 

Nj (s) N, (s) N3 (s) Nl (s) 0 
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equals the numerator polynomial exactly. But, this is simply a determinantal equation involving a 
matrix of polynomials and hence is solvable via the QZ-method and algorithm presented in the 
previous section. The validity of the fact that the determinant defined by equation 32 equals the 
numerator polynomial for the partial fraction expansion exactly defined by equation 31 can be 
easily verified by expanding the determinant in terms of the cofactors ofNj (s),N 2 (s), • ♦ • , (s). 

CONCLUSION 

An attempt has been made to summarize a large body of numerical error analysis literature in 
a language familiar to controls engineers. If your computing problems are of moderate-to large- 
order, beware. The algorithms that you now use with confidence for low-order systems may not be 
applicable. 
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APPENDIX A 


THE METHOD OF LEVERRIER 




THE METHOD OF LEVERRIER 


Use of good engineering judgment normally leads to a well-conditioned computing problem with 
good blocking characteristics. Consider the need to compute the matrix (s {] - B)'* , where B is 
assumed low order, equal to n. The method of Leverrier is given as follows (see Reference 14, 
page 166): 

Step 1. Write the adjoint of the matrix (S d -B) as: 

(sd-B)'^ = s"-‘ -C, s"-2 + • • • ±C„ , 

Step 2. Make use of the equations 

Co= 

c; = r,d-c, B 
C, = 7,d-C, B 
C„., = '(n-.t! B 

C„., B = d 


and 


trace [C.B] = (i + 1) j = o, i n-i 

to compute all coefficient matrices C^,, , . . . , ^ . Tliis can be done iteratively by first using 

Cq to compute 7 ^ , then Cj , 7 ^ , , et cetera. 


Step 3. Make use of the fact that all B are quasi-upper triangular, and let: 


be the n eigenvalues of B. The desired matrix inverse is thus given by: 


(sd-B)-* 


C„s"-‘ -C, s"-2 +' • • ±C^ j 
(s - s, ) (s - s^) * • • (s-s^) 


Again, all B are of low order and hence the power series form of polynomials is acceptable. 
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